UNILE-CBR-2001-3 
OUTP-01-33P 



SUSY QCD and High Energy Cosmic Rays 1. 
Fragmentation functions of SUSY QCD 

Claudio CorianoQ and Alon E. FaraggiQ 

1 Dipartimento di Fisica 

Universita 7 di Lecce 
I. N.F.N. Sezione di Lecce 
Via Arnesano, 73100 Lecce, Italy 

2 Theoretical Physics Department 
University of Oxford, Oxford, 0X1 3NP, United Kingdom 

Abstract 

The supersymmetric evolution of the fragmentation functions (or timelike 
evolution) within N = 1 QCD is discussed and predictions for the fragmenta- 
tion functions of the theory (into final protons) are given. We use a backward 
running of the supersymmetric DGLAP equations, using a method developed 
in previous works. We start from the usual QCD parameterizations at low 
energy and run the DGLAP back, up to an intermediate scale - assumed to 
be supersymmetric - where we switch-on super symmetry. From there on we 
assume the applicability of an N = 1 supersymmetric evolution (ESAP). We 
elaborate on possible application of these results to High Energy Cosmic Rays 
near the GZK cutoff. 
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1 Introduction 



The Standard Model of elementary particle physics successfully accounts for all the 
observed particle properties. Its remarkable success suggests the possibility that 
it correctly describes the particle properties up to energy scales which far exceed 
the currently accessible energy scales. Indeed, this possibility has been entertained 
in the explorations of Grand Unified, and Superstring, Theories. Additionally, the 
spectacular confirmation of the Standard Model in high energy colliders, as well as 
the proton longevity and suppression of left-handed neutrino masses, provide strong 
support for the grand desert scenario and unification. The question then is how are we 
going to find out whether this is indeed the path chosen by nature, as direct signatures 
at LHC/NLC will be able to probe the desert only up to a few TeV. Therefore we 
need new tools to achieve this. An example of such a probe that has been used with 
great success is that of the gauge coupling unification. The general methodology in 
this respect is to extract the gauge couplings from experiments that are performed at 
the accessible energy scales. Based on concrete assumptions in regard to the particle 
content in the desert, the couplings are then extrapolated to the high energy scale and 
the unification hypothesis is tested. In this manner the consistency of the specific 
assumptions in regard to the particle content in the desert with the hypothesis of 
the gauge coupling unification is subjected to an experimental test. The well known 
spectacular success of this methodology is in differentiating between gauge coupling 
unification in supersymmetric and non-supersymmetric Grand Unified Theories [|T|. 
The unification hypothesis is found not to be consistent with the low energy data, 
unless the minimal particle content is supersymmetric, or is modified in some other 
way. Our task then is to develop additional tools that can serve as similar useful 
probes of the desert. 

1.1 Motivations for this work 

In this paper we start applying this philosophy to another physical setting. 

A different experimental observation - which also points toward the validity of the 
big desert scenario - is the existence of Ultra High Energy Cosmic Rays above the 
Greisen-Zatsepin-Kuzmin cutoff. A plausible explanation for the observation of such 
cosmic rays is the existence of supermassive metastable states with mass of the order 
10 12 — 10 13 GeV, and a lifetime that exceeds 10 w y. The possibility that such states 
compose a substantial component of the dark matter as well as that they may explain 
the observation of UHECR has been discussed elsewhere. It has been further shown 
that realistic string models often give rise to states with precisely such properties 0. 

The existence of such super-massive, slowly decaying states, may therefore serve 
as a probe of the grand desert. A plausible explanation of the observed showers in 
the UHECR is that the original supermassive decaying particle decays into strongly 
interacting particles which subsequently fragment into hadrons ||. This hypothesis 
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therefore needs the relevant fragmentation functions at the energy scale of the origi- 
nally decaying particle. These functions, however, similarly to the case of the gauge 
couplings considered in much of the literature on unification, are extracted from ex- 
periments at the currently accessible (low) energy scales which are much below that 
of the decaying particle. It is obvious that this points toward the use of the renormal- 
ization group (RG) evolution. However, this evolution, as well as the composition of 
the produced shower, depends on the assumptions made in regard to the composition 
of the spectrum in the desert. We conclude that similarly to the gauge coupling uni- 
fication methodology, the evolution of the fragmentation functions by utilizing the 
renormalization group extrapolation, and the subsequent implementation in the anal- 
ysis of the UHECR hadronic showers, can serve as a tool that differentiates between 
different assumptions on the very high energy spectrum. 

With these motivations in mind, here we develop the instruments to evolve the 
fragmentation functions in supersymmetric QCD up to a high energy scale. In par- 
ticular, in this paper we present for the first time the structure of the supersymmetric 
DGLAP (Exact Supersymmetric DGLAP or ESAP) equations in the timelike region. 
We solve them numerically for all flavours and introduce various combinations of 
non-singlet matrix equations to reach this result. We then study the impact of the 
new evolution assuming a common scale for SUSY breaking and restoration, param- 
eterized by the mass of the superpartners, here assumed to be mass degenerate. 

1.2 Step Approximations 

Given the long stretch it takes to proceed with the analysis of the SUSY evolution, a 
topic of numerical complexity in its own, the applications to cosmic rays of our work 
will be presented in a companion paper that will follow shortly. We also have decided 
to address issues related to the fragmentation region of N — 1 QCD starting from 
the "low energy" end, since nothing is known of these functions at any scale. The 
term "low energy" (the results we discuss here are obtained in the 10 3 — 10 8 GeV 
energy range ) - from the point of view of collider phenomenology - is a misnomer, 
but it is not in the context of high energy cosmic ray physics. There is also another 
more direct reason for limiting our analysis to this range. We have discovered some 
features of the RG evolution which require a special care at such large energies. We 
find -in our numerical studies- the onset of an instability in the evolution in the 
singlet sector of the fragmentation which requires an independent investigation. It 
is not clear to us, at this point, whether this instability is an intrinsic limitation of 
the algorithm used in the numerical implementation or if the perturbative expansion 
needs a resummation. At such large scales this last option remains open. We should 
also mention that we work in "x" space, and this might be a limitation. We hope to 
return to this issue in the near future. 

The analysis of fragmentation functions of SUSY QCD is a topic of remarkable 
phenomenological interest both for collider phenomenology and in the astroparticle 



3 



context. While - recently - a detailed analysis of the supersymmetric evolution of 
the parton distributions of SUSY QCD has been presented ||, |], the study of the 
evolution of the supersymmetric fragmentation functions is still missing. Of partic- 
ular interest is the study of a combined QCD/SQCD evolution with intermediate 
regions in the evolution characterized by partial supersymmetry (SAP) or exact su- 
persymmetry (ESAP) ||. It is well known that a matching between these regions is 
possible with specific boundary conditions. These boundary conditions should come 
from some deeper understanding of the way in which supersymmetry is broken and 
restored as we run into the different stages of the renormalization group evolution. 
In a first approximation, however, it is possible to assume that the regular QCD dis- 
tributions (or fragmentation) functions are continuous at each intermediate region, 
thereby neglecting threshold effects which cannot be inferred from first principles. As 
we cross any supersymmetric region, starting, for instance, from the end of the reg- 
ular QCD evolution, it is possible to assume that supersymmetric distributions and 
fragmentation functions are generated radiatively. We recall that a similar modeling 
-which neglects possible effects of higher twist at the thresholds- is well spread in the 
case of ordinary QCD. 

2 Initial/Final State Scaling Violations 
2.1 The Spacelike Evolution 

The first analysis of the spacelike evolution was carried out in ref. [[J, where simple 
models -obtained from the analysis of the first 2 moments - of the supersymmetric 
parton distributions were also presented. A complete strategy to solve these equations 
and the strategy to generate supersymmetric scaling violations in N=l QCD was put 
forward in Refs. [f|, |5[. 

In N = 1 QCD gluons have partners called gluinos (here denoted by A) and left- 
and right-handed quarks have complex scalar partners (squarks) which we denote as 
q"i and q~ R with q = q L + q R (for left-handed and right-handed squarks respectively). 

The interaction between the elementary fields are described by the SU(3) color 
gauge invariant and supersymmetric Lagrangian 

+iqi pq { + D^q R D»q R + D^q L D»q L + igy/2 (\ a R q\ L T a q Li + \ a L q\ R T a q Rl - h. c.) 
-\d 2 (q[i Ta QLi - q Ri T a q R i) 2 + mass terms, 

(1) 

where a runs over the adjoint of the color group and i denotes the flavours. 
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qm are the supersymmetric partners of right-handed quarks qui and qn are those 
of the left-handed quarks qn. 

We remind here that in actual phenomenological applications, it is useful to ex- 
tract the light-cone dynamics -or parton model picture- of supersymmetric collisions 
by invoking a factorization of the cross section into hard and soft contributions. Frag- 
mentation functions appear in this context. An integral part of this analysis is the 
use of renormalization group equations -SUSY DGLAP- which resum the logarithmic 
scaling violations and evolve distributions functions and fragmentation functions to 
the appropriate factorization scale. 

In our approach R-parity is conserved. However, we neglect contribution from 
fragmentation into stable supersymmetric states which are subleading relative to 
fragmentation into strongly interacting states. Similarly, we do not consider the 
possibility of R-parity violating operators. Limits on such operators [|7| would suggest 
that they cannot produce sizable effect on the fragmentation function. However, more 
detailed studies may be required to ascertain this general expectation. 

These studies require the matrix of the anomalous dimensions (for all the mo- 
ments), or of the corresponding Altarelli-Parisi (DGLAP) kernels. We recall that 
the leading order anomalous dimensions are known || || both for a partial SUSY 
evolution and for the exact one. 

2.2 The Timelike Evolution 

In this work we focus our attention on the evolution of fragmentation functions in 
the context of exact supersymmetry, with coupled gluons and squarks. 

Fig. 1 summarizes the basic strategy of our work. The figure at the top depicts 
a regular evolution (a one-phase evolution), while the one at the bottom illustrates 
a mixed QCD/SQCD evolution (a two-phase evolution). The first stage is super- 
symmetric, the second one is regular. In the first stage all the partons (SUSY and 
not SUSY) are massless, in the second one only quark and gluons survive. At the 
end of the first stage the supersymmetric partners are extinct, but their presence 
at the higher scales is saved in the boundary conditions for the quarks and gluons 
fragmentation functions, when the second stage of the evolution starts. 

Within the approximation implied by the use of a pure SQCD evolution, no 
missing E t events are allowed along the development of the supersymmetric cascade, 
since electroweak and flavour mixing effects are not included. Their inclusion is still 
an open chapter even in (QCD) Monte Carlo event generators for the final state. 

The ladder structure of the supersymmetric showers in the leading logarithmic 
approximation can be easily pictured in terms of symmetric doubling of lines of 
regular QCD ladders in all the possible allowed cases. Showers can be equally initiated 
by q q pairs, gluon pairs, gluino pairs or squark pairs. The supersymmetric transition 
(i.e. those involving supersymmetric partners) stop (see Fig. 1 (bottom)) once the 
supersymmetry breaking scale M\ = Mq is reached. 
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Similarly to the QCD case, in the case of exact N — 1 supersymmetry we de- 
fine singlet and non-singlet fragmentation functions Dj(x,Q 2 ). They describe the 
amplitude for a parton of type / to fragment into a hadron h as a function of the 
Bjorken variable x (fractional energy of the fragment) and initial energy Q. Their 
operatorial definition is similar to that of ordinary parton distributions. The evolu- 
tion equations for these functions are similar to the standard DGLAP equations, but 
with transposed kernel matrix (P — > P T ). In leading order the analytic continuation 
(from the spacelike or DIS evolution to the timelike evolution) is straightforward, 
while some complications appear to next-to-leading order. They involve a breaking 
of the Drell-Yan-Levy relation, which appears to be violated at parton level |J. Since 
our analysis is in the leading logarithmic approximation, we will not consider these 
aspects in our work any further. 

The equations for the timelike evolution are given by 



d 



dhg(Q 



d 

d\og(Q 2 ) 
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D h x (x,Q 2 ) 
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d\og(Q 



d 



D h ,A^Q 2 ) 



dlog(Q 2 ) mL 
- D h ~ (x O 2 ) 
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2^ 
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2tt 
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p 99 ® D h g + p Xg ® D h x +p qg ®Yi { D i + D i 
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n s _ _ \ 

+Pqg ® (fti + + + QiR) 
t=l / 

P gX ® Pj + P\\ ® P A + P 9 a ® £ K + D i 



(2) 



"1 

i=l 



(3) 



1 



P A9 ® p£ + p^ g 



An, 



9 ' An 



+ p m ® 



4n 



An, 



-Pw ® D\ 



In the following we will use the short-hand notation 

pJ>,g 2 ) = p£ L (*,g 2 ) + p^(x,g 2 



p~ +p 

9iL 



/i 



^ ® + — Pxa ®D h ,+ -P ad <g> P^ 



(4) 



(5) 



(6) 



(7) 



to denote the fragmentation functions of squarks of flavour i at a fractional energy 
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x and momentum Q. It is also convenient to separate the equations, as usual, into 
singlet and non singlet sectors using the definitions 



D h qv (x,Q 2 ) = eK(^v^,(^ : 



D h 

QNS 



1=1 



D q+ (x,Q 2 ) = J2{D k qi (^Q 2 ) + Di(^Q 2 )) 



i=i 



i=l 



The non singlet equations are 



Q 2 ^Dq v (x,Q 2 ) 
d 



Q W 2 ^ M) = 



«(g 2 

a(Q 2 
2tt 



(8) 



(9) 



and the singlet matrix equations, which mix qy and qy with the gluons and the 
gluinos 
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r qq J 



® 



D*(x,Q 2 ) 
d\(x,Q 2 ) 
D h M (x,Q 2 ) 

D h gm (x,Q 2 ) 



(10) 



where "T" indicates the matrix transposed. To solve for all the flavours, it 
nient to introduce the linear combinations 



is conve- 



Xi (x,Q 2 ) = D\ +) 



1 

n .f 
1 



D h qM 



Ux,Q 2 ) = D% ) -—D h 5{ 



(+) 



(11) 



and the additional singlet equations 
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and 



dQ 2 ^ > 2tt v qq *r qq q 



(-) 



(12) 



(13) 



The general flavour decomposition is obtained by solving the singlet equations for 
D h (+) and D\ +) , then solving the non-singlet equations for -D h ( _ } and D\ +) and for 

and D% . The fragmentation functions of the various flavours are extracted using 
the linear combinations 



(14) 

and 

(15) 

to identify the various flavour components. 

There are simple ways to calculate the kernel of the SUSY DGLAP evolution by 
a simple extension of the usual methods. The changes are primarily due to color 
factors. There are also some basic supersymmetric relations which have to be sat- 
isfied ||. They are generally broken in the case of decoupling. We recall that the 
supersymmetric version of the (3 function is given at two-loop level by 
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(31 = -{llC A -2n f -2n x ) 

$ = ^(3AC 2 A -10C A n f -10C A n x -6C F n f -6C x n x ) (16) 

where n/ is the number of flavours and for Majorana gluinos n x = C A - 

The ordinary running of the coupling is replaced by its supersymmetric running 

a S {Ql) = 2 1 / _ /?f lnln(Q 2 /A 2 ) 1 A 

2vr /3 5 ln(Q 2 /A 2 ) ^ $ + l ln 2 (g 2 /A 2 ) 7 ' 

The kernels are modified both in their coupling (a —>■ a s ) and in their internal 
structure (Casimirs, color factors, etc.) when moving from the QCD case to the 
SQCD case. In order to illustrate the approach that we are going to follow in our 
analysis of the fragmentation functions within a mixed SQCD/QCD evolution, we 
recall the strategy employed in Refs. H to generate the ordinary supersymmetric 
distribution functions. 

In Ref. |J for scaling violations affecting the initial state were introduced 3 regions: 

• 1) the QCD region, described by ordinary QCD (Altarelli Parisi or AP); 

• 2) an intermediate supersymmetric region with coupled gluinos and decoupled 
squarks (partially supersymmetric AP or SAP); 

• 3) the N=l region (exact supersymmetric AP or ESAP). 

The basic idea was to generate these distributions radiatively as is usually done in 
QCD for the gluons, for instance, using the fact that the matrix of the anomalous 
dimensions is not-diagonal. The most general sequence of evolutions (denoted AP- 
SAP-ESAP in 0) is described by the arrays (Qi, Qf)Ap — (Qi, Qf)sAp — (Qi, Q/)esap, 
with Qf,AP = Qi,sap = Tfi2\ and Qf,sAP = Qi,ESAP = friq- In this work we limit our 
analysis to a simpler AP-ESAP evolution and we run the DGLAP equations back 
starting from known QCD fragmentation functions - for which various sets are avail- 
able in the literature |TT[ - up to an intermediate supersymmetric scale. From this 
point on (up in energy) we switch on the ESAP evolution of fragmentation functions. 
As we run the equations upward, supersymmetric fragmentation functions are gener- 
ated. The boundary values of the low energy (QCD) functions at the supersymmetry 
scale set the initial condition for the (backward) supersymmetric running up to 
the final scale. In simple terms: we reach the mountain from the valley, and cross a 
fence along the way. 

In general, if we split the intermediate SUSY scale into 2 sectors m x « rriq 
(ESAP-SAP-AP evolution) in this (general) case the solution is built by sewing the 
three regions as 
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D h qv (x,Q 2 ) 
Dl(x,Q 2 ) 



+ f 2X d\og Q 2 P^ P NS (x,a(Q 2 )) 
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(18) 



in the non singlet and 
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J mi. 



D h G (x,Q 2 ) 
D h x (x,Q 2 ) 
D h q+ {x,Q 2 ] 




+ 



f Q * f d\ogQ 2 P» s ESAP (x,a ES (Q 2 ))® 

J mi. 



D G (x,Q}) 
D h x (x,Q 2 ) 

D h q+ (x,Q 2 ) 
D h 2q+ (x,Q 2 ) 



(19) 



for the singlet solution, where rri2 q = 2m q . The zero entries in the arrays for some 
of the distributions are due to the boundary conditions, since all the supersymmetric 
partners are generated, in this model, by the evolution. 

The general structure of the algorithms that solves these equations is summarized 
below. We start from the non singlet sector and then proceed to the singlet. From now 
on, to simplify our notation, we will omit at times the index "/i" from the fragmen- 
tation functions when obvious. We define D q Ns(x ) Q 2 ) = (D qv (x,Q 2 ), D qv (x,Q 2 )) 

and DA n (x) = (DA q v ,DA^ and introduce the ansatz 



n 



D qNS (x,Q 2 ) = J2 

n=0 




(20) 



where n is an integer at which we stop the iteration, which usually ranges up to 20. 
The first coefficient of the recursion is determined by the initial condition 



DA (x) = U(x)®D" s (x,Ql), 



(21) 
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where 



U x {x) = 
U 12 (x) = 



5(1 -x) 



5(1 -x) 

5(1 -x) 



(22) 



The recursion relations are given by 



DA 



n+l 



Po 



T,NS 
AP 



DA n (x) 



(23) 



by 



The solution in the first (DGLAP) region at the first matching scale m 2 \ is given 



D q Ns(x,m 2 \) 



"0 f)AAP( \ 
n=0 



'a(m 2 \ 2 ) 
, a(Ql) , 



(24) 



valid for m\ < rriq. As we will discuss in the next sections, these two scales will 
be taken to be degenerate and the numerical treatment simplifies considerably. We 
also mention that we don't find any significant difference, from our analysis, if this 
degeneracy between the gluino and squark masses is lifted. However, just for illus- 
trative purposes, we keep the two scales separated in the formal analysis presented 
below in order to show how the matching is performed in a more general case. Mass 
differences within a range of 100-200 GeV play a rather modest role in the analysis 
and are hardly distinguishable. 

At the second stage the (partial) supersymmetric coefficients are given by (S is a 
short form of SAP) 



DA s ' NS (x) = U(x)®D q (x,m 2X 2) 

DA^ s (x) = *Pj> NS (x)®DA^(x). 
Po 



(25) 



We construct the boundary condition for the next stage of the evolution using the 
intermediate solution 

evaluated at the next threshold m 2 \. 

n^^l^j^m. (27) 



n=0 



n\ \a(m 2 \ 2 ) 
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The final solution is constructed using the recursion relations 



DA ES (x) = U(x)®q(x,mq) 



(28) 



DA E n s +1 {x) = - w - s P^{x)®DA E n s {x) 



The final solution is written as 



D q Ns{x,Q) = D q Ns (x, Q Q ) + — - — log — — — - +^ — log' 



n\ \ a(Ql) ) n! V« 5 ( m 2A 2 ), 



n =i n! \a £5 (m^ 2 ) / 

(29) 

As we have already mentioned above, in the analysis presented below the two 
SUSY scales will be collapsed into one (m 2 \). 



3 QCD Evolution of the Fragmentation Functions 

As we have already mentioned, the timelike and the spacelike evolution, in leading 
order, are essentially the same. The singlet kernels are just the transposed of the 
kernels describing the evolution of the parton distributions. The ordinary QCD 
kernels, in leading order, are given by 



(0) ( l+x 2 \ 
p(0)( r ) - p(°) 

r qq \ X ) — r qq,NS 

PW(x) = n f (x 2 + (1 - xf) 
P«?(*) = C F l + { \- x)2 

pW(x) = 2N c (^^ + ^-2 + x(l-x)j+^8(l-x) (30) 

with p = 11/3 Ca — 4/3 TrU/ being the first coefficient of the QCD /5-function 
and Tr = 1/2. rif is the numbers of flavours. The equations for the fragmentation 
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functions in QCD are given by 



d\og(Q 2 ) q '" ' 2tt V 2n/ 



(31) 



and, as usual, can be decomposed into a non-singlet and a singlet sector 



Q 2 ^- 2 D h q Ux,Q 2 ) = ^-P qq , NS (x,a(Q 2 ))^D h q Ux,Q 2 ) (32) 



for the non-singlet distributions and 

\ / P.., P„ \ 



d 








dlog(Q 2 ) 




)■(?: 





D h qi+) (x,Q' 
D h (x,Q 2 ) 



(33) 



for the singlet sector, where 

M_ ) (* I Q a ) = l£-.D*. (34) 

A fast strategy to solve these equations, as discussed in ||, where a complete leading 
order evolution has been implemented, is to solve the recursion relations numerically. 
In the timelike case -that we are considering- these relations are obtained as linear 
combinations of the spacelike ones. We get 

n^"' a Cf f 1 dy yDAj \y) - xDJ^ \x) C f > n .,(-), 

DA^ +1 = -4—/ 4— log{l -x)DA% [x) 

Pq Jx y y-x (J 



+ 2^f f d JL (1 + z)DA i\ y) .^ D ^\ x) 

Po Jx y po 



(Jo Jx y y-x (J 

Po J a: y Po 

(35) 

A similar expansion is set up for the non-singlet variable Xi = Qi ~ V n / 1^ 

dazu = -i c ± f «L?E^m - -^m _ fjz , og(1 _ X)DA${X) 

(J Jx y y-x (Jo 
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+2 Cf f 1 ^V {l + z) DA x i{y) _ 3?fDA%(x). 
Po Jx y po 

(36) 

Similarly, the singlet equations generate recursion relations of the form 



DAi l\ (x) = /' ^ y^tM^tM _ 4 ^ log(1 _ x)DA f\ x) 

Po Jx y y-x Po 

+2 £l [ ld JL (1 + z)DA f> (y) _ 3 £f DA f> (x) 

Po Jx y Po 

Po Jx y y-x fa 

-2% f ^ (f + (1 - zf) DAf\y) - DA 9 n (x) 
Po Jx y 

-4^f ^(±-2 + *(l-*)W(y). (37) 



The coefficients of the various flavours both for and D% are then obtained 

yi Hi 

from the relations 



^ = l{ DA n+^- f DA i' ) ) 

DAI = \(dA% - ^-DA^ (38) 

4 The ESAP (N = 1 QCD) Fragmentation 

Moving to N — 1 QCD, we introduce recursion relations for appropriate linear com- 
binations of non-singlet fragmentation functions 

fa Jx y y-x po , 

-2%LDA?(x) - 2^f f %DA%(y) + 2^f f ^ (1 + z) DAf(y) 
Po Po Jx y Po Jx y 

(39) 
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-2^DA^{x) - 2^f f ^-zDA%(y) + 2^f f ^ (1 + z) DA%(y) 
Po Po Jx y p Jx y 

(40) 

and two similar non-singlet equations for D h ,_) and D h ( s 

-2^W-'(,) - 2^ /' *^-(v) + 4 /' * (1 + *) D^-'tv) 
Po Po A y Po Jai y 

(41) 



DA& = -4^ r ^ ggjtSO - ggjTw _ ^ lQg(1 _ x)d4%) 
p J x y y — x Po 

-2-fDA% (x)-2-f -^zDA% (y) + 2-f -H- (1 + z) DA% (y). 
Po Po Jx y Po Jx y 

(42) 

Analogous equations are satisfied by the combinations D^-) and D ! } ( _ ) by replac- 
ing in (H) and §|) L>V) D g(-) and £>-(-) D ^)- 
Moving to the singlet sector we obtain 



DAUM = _ 4 §1 /' ^ y£mH - ^M£> _ 4 £l log(l - X )DAI( X ) - DA i(x) 
Po Jx y y-x Po 

+ 2^d /' * (1 + ,) DA-(y) - f d JL (1 + , _ 2 ) M » W 
Po Ay po Jx y \z / 

+2 §i /■ ^ + (i - ,)•) !,>!.(,,) -2^ r^c^+d- 

Po Jx y ^ ' Po Jx y v 7 

_ 2 Ca ^ ^ + (1 _ ^ DA , (y) _^L( ld JL z{l _ z) DAUy) 
Po Jx y po Jx y 

(43) 

Po Jx y \z J Po Jx y y-x 
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-£± log(l - x)DA&x) + 2§i f ^ (1 + z) DAl{y) - DA*i 
Po Po Jx y 

Po Jx y Po 



x) 



_ 2 r}± f l djj (1 _ z) DAl(y) _ ^ ^ zDAi[y) (44) 
Po Jx y Po Jx y 



DAt + \{x) 



-2-f / + * 

Po •/* y \z 

(+) , 



2) DAi(y)-2^f ( ld f{l-z)DAl(y) 

/ Po J 2 y 



Cf ^dyyDAt\y)-xDAi 
Po 



,(+), 



2/ 

(/.;/ 



y- x 

7 (+) , 



2/ 

^logfl-^if'W 

Po 



+2^ / ^ (1 + z) DAf (y) - 2^ZMf ) (x) - 2^ / ^zl^ (y) 
Po •/* y Po Po Jx y 

(45) 



dy 



Po Jx y \z J Po Jx y ^ 1 

_ 4 CV ^dy yDAi{y)-xDAj{x) | ^ [ X *l DA *(y) 
Po Jx y y-x Po Jx y 

-4^r Ml " *)^f } (v) " ) (x) (46) 

Po Po 

where z = x/y. 

5 Numerical Results 

As an illustration of the procedure we adopt in our studies, let us consider the decay 
of a hypothetical massive state into supersymmetric partons. The decay can pro- 
ceed, for instance, through a regular qq channel and a shower is developed starting 
from the quark pair. The N — 1 DGLAP equation describes in the leading loga- 
rithmic approximation the evolution of the shower which accompanies the pair, and 
we are interested in studying the impact of the supersymmetry breaking scale (m\) 
on the fragmentation. In our runs we have chosen the initial set of Ref. [O], with 
parameterizations that can be found in the appendix. 

In our analysis we focus on the proton fragmentation functions and, following 



Ref. [n], we introduce the scaling variable 



In y' (47) 
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In LO we have = 88 MeV JTTJ and use three different values for /x , namely O] 



y/2 GeV, if a = u,d,s,g 

//q = <{ m(r? c ) = 2.9788 GeV, if a = c . (48) 

m(T) = 9.46037 GeV, if a = 6 

This leads to three different definitions of s. For definiteness, we use the symbol s c 
for charm and for bottom along with s for the residual partons. We parameterize 
the f.f.'s as 

D(x, fi 2 ) = Nx a (l - xf (l + ^ (49) 

and express the coefficients N, a, /3, and 7 as polynomials in s, s c , and s&. For 
s = s c = Sf, = 0, the parameterizations agree with Eq. (2) of Ref. [11| in combina- 
tion with the appropriate entries in Table 2 of that paper. The charm and bottom 
parameterizations must be put to zero by hand for s c < and s& < 0, respectively. 

Typical fragmentation functions in QCD involve final states with p, p, 7r ± , ir° and 
kaons k and we refer to the original literature for a list of all the fragmentation 
sets. We have chosen an initial evolution scale of 10 GeV and varied both the mass 
of the SUSY partners (we assume for simplicity that these are all degenerate) and 
the final evolution scale (Qf) in our backward evolution. At the end of the evolution, 
once the supersymmetric fragmentation functions are built, we interpret Qf to be 
the energy available for the decay of the massive state. We perform two types of 
investigations: 1) we analize the dependence of the evolution on the gluino mass m\, 
assumed to be degenerate with the quark mass; 2) we investigate the variation in size 
of the fragmentation functions in terms of the scale Qf. Since R-parity is conserved, 
the fragmentation scale is m 2 A, two times the gluino mass. Since our concern is in 
establishing the impact of supersymmetric evolution and compare it to standard QCD 
evolution across large evolution intervals, we plot the initial fragmentation functions, 
the regularly evolved QCD functions and the SQCD/QCD evolved ones. The latter 
two are originated from the same low energy form of Ref. [TT| . 

It is possible to include in the evolution of the fragmentation functions also differ- 
ent thresholds associated with more complex spectra in which the SUSY partners are 
not degenerate. These effects are negligible. Also, threshold enhancements may re- 
quire a more accurate treatment and will be discussed elsewhere. On general grounds, 
however, we don't expect them to play any important role, especially since we are 
interested in very extended renormalization group runnings. 

Let us now come to a description of our results. Fig. 2 shows the initial condition 
for the up quark fragmentation function into protons, and we have chosen the initial 
evolution scale (Qo) to be 10 GeV, as it has been extrapolated from collider data in 



Ref. [pT| . The evolution of this function follows QCD from this lowest scale up to 
a scale of 200 GeV (rri2\ = Qi), above which we use the full N = 1 evolution. In 
this figure we have chosen a lower initial fragmentation scale (Qf) of 1 TeV. As one 
can observe from the two plots for the up and down quarks, in general, the small-x 
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( diffract ive) region gets slightly enhanced when SUSY effects are taken into account, 
although the changes are quite small. The reason for this behaviour is to be found in 
the fact that the highest scale is not large enough to allow a discrimination of SUSY 
effects from the non supersymmetric ones, since not enough room is available for the 
super symmetric evolution. 

In Fig. 3 an analogous behaviour is found for the fragmentation functions of charm 
and strange quarks. Again, small-x enhancements are seen, but regular and SUSY 
evolution are hardly distinguishable. The situation appears to be quite different for 
the gluon fragmentation functions (f.f's) (Fig. 4). The regular and the SQCD evolved 
f.f.'s differ noticeably in the diffractive region, with a gluon fragmenting much faster 
when SUSY is in place, compared to the non supersymmetric case. This may be 
related to the wider phase space available for gluons to decay in the presence of 
supersymmetric channels. Qf is not too large, in this case, and equal to 1 TeV, 
which renders the plot particularly interesting from the experimental viewpoint. 

However, it should be kept in mind that even if the evolution predicts a larger 
probability for the fragmentation of a given parton into a specific hadron -compared 
to QCD-, this change does not automatically translates into a noticeable effect on the 
final multiplicity of that hadron in the final state. In fact, as we are going to show in 
the next section, at least in the case of the multiplicities, SUSY effects remain small 
up to pretty large fragmentation scales. 

As we raise the final evolution scale we start seeing some interesting features of 
the supersymmetric distributions. For instance, this is illustrated in Fig. 5 where we 
show the squark f.f.'s for all the flavours and the one of the gluino for comparison at 
a nominal fragmenting initial scale of 1 TeV. The scalar charm distribution appear 
to grow slightly faster than the remaining scalar ones, but the gluino f.f. is still the 
fastest growing at small-x values, which indicates a larger proability for gluino to 
fragment compared to scalar quarks. 

Overall, these probabilities, even at small-x values, remain small by a factor of 100 
compared to the ordinary QCD fragmentation functions. This is not a big surprise, 
given the low scale used in this specific example for Qf. 

In Figs. 6 and 7 we come to the second part of our analysis, and we start in- 
vestigating the dependence of the fragmentation functions on a varying initial scale 
Qf. We have selected the up quark and the gluon in these figures. The impact of a 
varying scale on these functions appear to be rather small in the quark case, while is 
more pronounced for the gluon case. Supersymmetry is broken and restored, in these 
two examples, at a scale of 200 GeV. 

An important point which deserves thorough attention is the study of the de- 
pendence of the fragmentation functions on the SUSY scale m\, which is shown in 
Fig. 8. We have varied this scale in a considerable range (from 200 GeV up to 1 TeV) 
and kept the initial fragmentation scale Qf fixed at 100 TeV. The changes for the 
up quark fragmentation functions are all within a 5 % interval. A similar result is 
presented in Fig. 9, where we illustrate in a same plot the dependence on the SUSY 
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scale m\ of the fragmentation functions of both the squark up and the up quark. 

The shape of the gluino density as a function of the the initial fragmentation 
scale is shown in Fig. 10. The growth appears to be rather pronounced at small-x. 
Increasing the fragmenting scale Qf, fragmentation occurs at larger x- values. 

Fig. 11 illustrates the behaviour of f.f's of all types at a large initial fragmentation 
scale. In order to illustrate the different behaviours that the quark/squark sector has 
compared to the gluon/gluino sector, we have shown the f.f.'s for the b quark and the b 
squark. Although all the functions get a small-x enhancement or are more supported 
in this x-region, the largest fragmentation appears to come from the gluon-gluino 
sector. The gluon fragmentation function is still the largest and fastest growing at 
small-x, due to the small-x behaviour of the gluon anomalous dimensions, while the 
fragmentation functions for the bottom quark and for its superpartner are rather 
small. 

Finally Fig. 12 has been included in order to clarify some issues concerning the 
possible impact of our results on the multiplicities of the final state hadrons due to the 
supersymmetric evolution. While a more detailed analysis of this and of other similar 
issues will be presented elswewhere, here we would like to make some comments in 
this directions. 

For this purpose, consider the decay of the metastable state into qq pairs (for 
QCD) and into a qq or a qq pairs (for SQCD ) mediated by a heavy photon (7*). In 
one case we follow a standard DGLAP evolution, in the second case we analize the 
probability for the supersymmetric partons to end up into protons using the SUSY 
evolution. In Fig. 12 we focus on the following observables, characterizing the final 
state 

R h Q c D {Q 2 )=T,^j zh dz(D h qi (z,Q")+D^z,Q 2 )) (50) 

and 

Rsqcd(Q 2 ) = E^ 2 / h dz {D h qi {z, Q 2 ) + D\{z, Q 2 ) + £>J.(z, Q 2 ) + D^z, Q 2 )) 

(51) 

with = rrih/(Q/2) being the minimum fractional energy required for the frag- 
mentation to take place, are the charges of the n/ quark and squark flavours, 

We recall that since Dp arton (x,Q 2 ) is the probability for a given parton to frag- 
ment into a final state hadron h at a given fractional energy x, given an original 
fragmentation scale Q, the two observables above describe the total probability for 
producing a hadron h, through quark channels and their corresponding supersym- 
metric counterparts. 

We study the numerical values of Rqcd(Q) and Rsqcd(Q) as a function of Q 
(in GeV), the initial fragmenting scale of partons into protons. We show in fig. 12 
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the variation of Rqcd and Rsqcd with the initial energy scale Q in a given interval 
(1 TeV - 4 TeV). The fragmentation probabilities, in this specific case, differ by 
approximately 2 %. The opening of supersymmetric channel, in our results, tends to 
slightly lower the probabilities for fragmentation -compared to standard QCD-. Here 
we have chosen a nominal gluino mass of 100 GeV. 

6 Conclusions and Perspectives 

We have presented a first detailed study of the fragmentation functions of SUSY 
QCD in the leading logarithmic approximation and presented, for the first time, 
results for the f.f.'s of all the partons into protons within a radiatively generated 
model. Our analysis here has been quite conservative and focused on an intermediate 
energy region where small-x effects, linked to resummation and Regge behaviour, 
which might invalidate a simple leading logarithmic (in Q 2 ) evolution, may still be 
under control. Within the obvious approximations involved in our study, we can still 
envision that some of the features present in our results may become more pronounced 
once we increase the original energy scale at which fragmentation occurs. However, 
some variants, compared to the ordinary QCD scenario, are noticed. 

An obvious question is why should we consider inclusion of supersymmetry in the 
fragmentation function in the first place. After all the fragmentation functions are 
measured at the GeV scale and current experiments requires their extrapolation to 
energy scales which are still below the SUSY threshold. In this respect our analysis 
was motivated from the recent exciting results from Ultra High Energy Cosmic Rays 
observations, which may indicate the existence of a new hadronic scale, of the order 
of 10 11 GeV, or so. If the UHECR results are substantiated further in the forthcoming 
experiments, as we elaborated in the introduction, it precisely necessitates the type 
of analysis that we pursued in this paper. It opens up the need to extrapolate the 
QCD parameters to the new scale, and the inclusion of any new degree of freedom 
which is crossed in the extrapolation. As we discussed here this also opens up a 
plethora of new issues and challenges for the QCD extrapolation and analysis. 

Although these variations are small at the energy probed in our numerical analysis, 
other sources of enhancement could set in at higher scales. First of all, at a higher 
fragmenting scale, we expect more remarked differences in the overall distributions 
of the total fragmenting multiplicities (above the 2 % level observed in Fig. 12), but 
we also expect that various (and less inclusive) correlations (such as energy-energy 
correlations and event shape variables) may affect more substantially the structure of 
a supersymmetric final state with respect to standard QCD. As we approach the GZK 
cutoff, these features may become much more pronounced, although the result has 
clearly to do with the distributions and coupling of the supersymmetric constituents 
in the metastable state that undergoes decay. 

The natural question to ask is then: what are the distributions of supersymmetric 
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partons in a decaying metastable state of very large mass prior to fragmentation? At 
such large scales, so far from the 1 TeV scale preferred by many standard supersym- 
metric models, the Renormalization Group strategy can be easily embraced to its 
fullest extent with sizeable consequences on the low energy end. On this and other 
related issues we hope to return in more detail in the near future. 
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7 Appendix 1. Input Functions 



For convenience we have included below the list of f.f.'s taken from Ref. [IT| that we 
have used 

Vlvi^ ,,2\. 



Z^(*,/i 2 ) = 2DfW 



D p J p (x,n 2 ): 



x, n 



N = 0.40211 - 0.21633s- 0.07045s 2 + 0.07831s 3 

a = -0.85973 + 0.13987s- 0.82412s 2 + 0.43114s 3 

(3 = 2.80160 + 0.78923s- 0.05344s 2 + 0.01460s 3 

7 = 0.05198s - 0.04623s 2 (52) 



N = 4.07885 - 2.97392s- 0.92973s 2 + 1.23517s 3 

a = -0.09735 + 0.25834s- 1.52246s 2 + 0.77060s 3 

P = 4.99191 + 1.14379s- 0.85320s 2 + 0.45607s 3 

7 = 0.07174s - 0.08321s 2 (53) 



N = 0.11061 - 0.07726s c + 0.05422s 2 - 0.03364s 3 

a = -1.54340 - 0.20804s c + 0.29038s 2 - 0.23662s 3 

P = 2.20681 + 0.62274s c + 0.29713s 2 - 0.21861s 3 

7 = 0.00831s c + 0.00065s 2 (54) 
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N = 40.0971 - 123.531s 6 + 128.666s 2 - 29.1808s 3 , 

a = 0.74249 - 1.29639s b - 3.65003s 2 + 3.05340s 3 , 

P = 12.3729 - 1.04932s 6 + 0.34662s 2 . - 1.34412s 3 , 

7 = -0.04290s 6 - 0.30359s 2 (55) 



N = 0.73953- 1.64519s + 1.01189s 2 -0.10175s 3 

a = -0.76986 - 3.58787s + 13.8025s 2 - 13.8902s 3 

= 7.69079 -2.84470s - 0.36719s 2 - 2.21825s 3 

7 = 1.26515s- 1.96117s 2 (56) 



8 Appendix 2. The Weights of the N=l Kernels 

We briefly recall the numerical strategy employed in this analysis. A more detailed 
description will be given elsewhere. We just mention that the radiative generation of 
super symmetric distributions requires special accuracy since these scaling violations 
grow up very slowly. We define P(x) = xP(x) and A(x) = xA(x). We also define 
the convolution product 

The integration interval in y at any fixed x-value is partitioned in an array of in- 
creasing points ordered from left to right (x , xi, x 2 , x n , x n+1 ) with x = x and 
x n+ i = 1 being the upper edge of the integration region. One constructs a rescaled 
array (x, x/x n , x/x 2 , x/xi, 1). We define Sj = x/xi, and s n+ i = x < s n < s n _i < 
...si < s = 1. We get 

At this point we introduce the linear interpolation 

A(y) =(l- V ~ X _ l ) A( Xi ) + V ~ X _ % A(x i+1 ) (59) 

and perform the integration on each subinterval with a change of variable y— > x/y 
and replace the integral J(x) with its discrete approximation Jn(x) to get 



J N (x) = A(x )-^— f ^P(y)(y - Sl ) 
1 - si Jsi y 
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+ J2M x i) — - — f djL P{y){y - s l+l ) 

„_1 Si Sj4_1 JsiJ-1 V 



i=l 
N 



s i s i+i Jsi+i y 

I ' : (ItJ 



i = i Si— i Si J Si y 



(60) 



with the condition A(xn + i) = 0. Introducing the coefficients W(x,x) and W(xi,x), 
the integral is cast in the form 



J N (x) = W(x, x)A(x) + E ^(*i, x)^(xi) 



i=i 



where 



1 


— Si 7si 




Si 




— s i+l 




Si 



Si-l — Sj ,/s; |/ 



— (y - Si+i) P(y) 
^(y-s^Piy). 



We recall that 



Jo 



dx- 



■ d!/ /W-/(i) 



(l-x)+ Jo l-y 



and that 



i ^ m ^[ 1 dy yf(y)--m +f{x)log{1 _ x) 



(l-x) 



x y y-x 



as can can be shown quite straightforwardly. 
We also introduce the expressions 



(61) 



(62) 
(63) 

(64) 



In (x) 

Jrii(x) 
Jnti(x) 
Jnt\ (x) 



1 



l-si 
1 



1 

Sj— 1 Sj 

1 



log(si) + log(l - Si) 



iog(i^ ±i ) + « + ik>g /r 1 ~ s * Sm 



log 



I - Si J 

' 1-Sj \ 

. 1 - s i-l J 



s 1 — Sj+i Sj 

1 ( Si \ , f l ~ s i-\ 
+ Sj_i log ( ) + Sj„i 



Si-l . 



1 - Si 



l-si 



log Si. 



, i = 2,3, ..JV 



(65) 
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Using the linear interpolation formula fl59|) we get the relation 



1 dyyAJy) — xAJx) , ,„ . . . . . . . T . . 

— — = -log{l-x)A n (x)+A n (x)In (x) 

x y y-x 

N 

+ A n (xi) (Jrii(x) - Jnti(x)) (66) 

8=1 

which has been used for a fast and accurate numerical implementation of the 
recursion relations. 
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Figure 2: Fragmentation function for the up (top graph) and down quark (bottom 
graph) xD^ p d (x, Q 2 ) at the lowest scale (input) Qq = 10 GeV, and their QCD (or 
regular) and SQCD/QCD evolutions with Qf =10 3 GeV. The SUSY fragmentation 
scale is chosen to be 200 GeV. 
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Figure 3: Fragmentation functions into protons of the charm and strange quarks 
xDP'P(x, Q 2 ) at the lowest scale (input) Qo = 10 GeV, and its evolved QCD (regular) 
and SQCD/QCD evolutions with Qf =10 3 GeV.The SUSY fragmentation scale is 
chosen to be 200 GeV. 




Figure 4: The gluon fragmentation function xD p ' p (x, Q 2 ) at the lowest scale (input) 
Qo = 10 GeV, and its evolved QCD (regular) and SQCD/QCD evolutions with 
Q f =10 3 GeV.The SUSY fragmentation scale is chosen to be 200 GeV. 
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Figure 5: The fragmentation functions of squarks and gluino at the lowest scale 
(input) Q = 10 GeV, with Q f = 10 3 GeV and SUSY scale 200 GeV. 
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Figure 6: The up quark fragmentation function xDP' p (x, Q 2 ) evolved with 
SQCD/QCD for varying final values of Qf. The SUSY fragmentation scale is chosen 
to be 200 GeV. 
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Figure 7: The gluon fragmentation function xD^' p (x, Q 2 ) evolved regularly and ac- 
cording to SQCD/QCD for varying final values of Qf. The SUSY fragmentation scale 
is chosen to be 200 GeV. 
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Figure 8: The up quark fragmentation function xD^' p (x, Q 2 ) evolved regularly and 
according to SQCD/QCD for varying final values of the SUSY scale m 2 \ and Qf = 10 5 
GeV. The SUSY fragmentation scale is chosen to be 200 GeV. 
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Figure 9: The quark and squark (up) fragmentation functions for a varying vrt\ and 
a fixed large final scale Qf = 10 5 GeV. 
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Figure 10: The gluino fragmentation function xD p,p (x, Q 2 ) evolved in the region 
200 — Qf GeV for different values of Qj The SUSY fragmentation scale is chosen to 
be 200 GeV. 
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Figure 11: The quark and squark (bottom) fragmentation functions together with 
those of gluon and gluino for m\ = 100 GeV and a fixed large final scale Qf = 10 5 
GeV. 
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Figure 12: Rqcd and Rsqcd for protons in the final state versus Q (in GeV). 
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